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ABSTRACT 

Problems  requiring  regression  analysis  of  censored  data 
arise  frequently  in  practice.  For  example,  in  accelerated 
testing  one  wishes  to  relate  stress  and  average  time  to  failure 
from  data  including  unfailed  units,  i.e.,  censored  observations. 

Maximum  likelihood  is  one  method  for  obtaining  the  desired 
estimates;  in  this  paper,  we  propose  an  alternative  approach. 

An  initial  least  squares  fit  is  obtained  treating  the  censored 
values  as  failures.  Then,  based  upon  this  initial  fit,  the 
expected  failure  time  for  each  censored  observation  is  estimated. 
These  estimates  are  then  used,  instead  of  the  censoring  times, 
to  obtain  a revised  least  squares  fit  and  new  expected  failure 
times  are  estimated  for  the  censored  values.  These  are  then 
used  in  a further  least  squares  fit.  The  procedure  is  iterated 
until  convergence  is  achieved.  This  method  is  simpler  to  imple- 
ment and  explain  to  non-statisticians  than  maximum  likelihood 
and  appears  to  have  good  statistical  and  convergence  properties. 

The  method  is  illustrated  by  an  example,  and  some  simulation 


results  are  described.  Variations  and  areas  for  further  study 


also  are  discussed. 


KEY  WORDS:  Censored  data,  regression  analysis,  maximum 
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INTRODUCTION 


Statistical  relationships  between  variables  must 
often  be  fitted  from  data  including  censored  observations. 

For  example,  in  evaluating  material  properties,  life  tests 
are  conducted  at  accelerated  stresses;  the  resulting  data 
are  used  to  relate  stress  and  average  time  to  failure 
assuming  a model,  such  as  the  Arrhenius  relationship  or  the 
inverse  power  law.  However,  at  low  stress  conditions  --  the 
region  of  greatest  interest  — some  test  units  remain 
unfailed  by  the  time  the  data  must  be  analyzed,  thus  resulting 
in  censored  observations.  Indeed,  a test  program  which  is 
designed  to  assure  that  all  units  fail  within  the  available 
testing  time  might  require  stresses  to  be  accelerated  to  such 
an  extreme  that  new  failure  modes  are  introduced  and  the 
assumed  model  is  no  longer  valid. 

Maximum  likelihood  is  one  method  for  obtaining  the 
desired  estimates  from  the  resulting  censored  data.  In 
practice,  however,  data  from  test  programs  involving  cen- 
soring are  frequently  analyzed,  not  by  a professional 
statistician,  but  by  an  engineer  or  scientist  who  has  had 
only  limited  training  in  statistics.  Such  persons  generally 
have  not  been  exposed  to  maximum  likelihood  methods  and 
may  not  have  easy  access  to  a computer  program  for  performing 
the  complex  required  calculations.  Instead,  we  have  found 
that  censored  data  are  frequently  analyzed  incorrectly  by 
standard  least  squares  techniques  for  complete  data  using 
one  of  the  following  two  expediencies  to  remove  the  censoring: 
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- Treating  the  censored  observations  as  if  they  were 
uncensored,  e.g.,  in  analyzing  product  life  data, 
assuming  the  censoring  times  to  be  failure  times. 

- Ignoring  the  censored  observations  altogether. 

Both  of  these  procedures  lead  to  biased  estimates  of 
the  regression  line.  In  particular,  they  result  on  the 
average  in  overly  conservative  estimates  of  product  life 
at  design  conditions;  the  first  procedure  incorrectly 
assumes  that  a failure  has  taken  place  when,  in 

fact,  it  has  not,  while  the  second  approach  ignores 
valid  data  on  unfailed  units  (often  at  the  conditions 
of  greatest  practical  interest) . These  two  incorrect 
procedures  also  produce  biased  estimates  of  the  variability 
around  the  fitted  line,  as  measured  by  the  residual  error 
variance.  Both  cases  underestimate  the  values  on  the 
average,  resulting  in  too  short  a confidence  interval 
on  the  average  life  at  a design  condition. 

It  is  therefore  desirable  to  have  methods  for 
analyzing  censored  data  which  are  relatively  simple  to 
explain  to  non-statisticians  and  to  apply  by  them,  in 
addition  to  having  good  statistical  properties.  One  such 
method,  based  upon  an  iterative  least  squares  approach, 
is  proposed  here.  To  apply  this  method  to  fit  standard 
relationships  for  censored  data  under  the  usual  assumptions, 
such  as  a normal  or  log  normal  distribution  for  the  random 
variation,  requires  only  a computer  routine  for  least 
squares  regression  analysis  and  a tabulation  or  calculation 
of  normal  distribution  ordinates  and  areas. 
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Most  of  the  subsequent  discussion  will  be  in  the 
context  of  product  life  analysis  and  thus  will  be  con- 
cerned with  censoring  to  the  right.  Also,  we  will  be 
principally  concerned  with  the  single  censoring  situation 
where  the  censoring  times  exceed  the  faj lure  times  at  each 
stress,  i.e.,  so-called  Type  I or  time-censoring.  As  we 
will  indicate  later,  the  basic  approach  also  applies  for 
more  general  situations. 

In  summary,  the  proposed  method  is:  At  Step  1 of 
Iteration  0 obtain  an  initial  least  squares  fit  from  the 
given  data  treating  the  censored  values  as  failures.  At 
Step  2 of  this  iteration  use  this  initial  fit  to  estimate 
the  expected  failure  time  for  each  unfailed  unit,  conditional 
upon  its  observed  censoring  time.  To  obtain  a revised  least 
squares  fit  at  Step  1 of  Iteration  1,  use  these  estimates  of 
expected  failure  times  for  the  unfailed  units  in  place  of 
the  censoring  times.  At  Step  2 of  Iteration  1 
estimate  from  this  new  fit  the  new  expected  running  times 
for  the  censored  observations.  At  Iteration  2 use  these  now 
estimates  in  a further  least  squares  fit.  Iterate  this  pro- 
cedure until  convergence  is  achieved. 

In  Section  2 of  this  paper,  relevant  past  results  are 
reviewed.  Section  3 provides  details  of  the  proposed 
iterative  least  squares  procedure  for  a simple  linear  model, 
normally  distributed  random  variation  and  censoring  to  the 
right.  This  method  is  illustrated  by  a numerical  example 
in  Section  4.  These  analyses  should  suggest  the  use  of  the 
method  for  more  general  situations.  In  Section  5 the  statis- 


tical  properties  of  the  proposed  procedure  are  evaluated 
by  a simulation  analysis  and  compared  with  those  of  the 
method  of  maximum  likelihood.  Generalizations  of  the 
results  and  areas  for  further  study  are  discussed  in 
Section  6 and  7 respectively,  and  some  concluding  remarks 
are  made  in  Section  8. 
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2 * RELATED  work  and  literature 

Hahn  and  Nelson  C 3 J describe  and  compare  three  general 
methods  for  analyzing  censored  life  data  to  estimate  rela- 
tionships between  stress  and  product  life: 

- Maximum  likelihood  estimation 

- Linear  unbiased  methods 

- Graphical  methods. 

The  method  of  maximum  likelihood  is  the  most  general  of 
these  approaches  and  is,  perhaps,  the  one  used  most  frequently 
by  statisticians.  However,  it  requires  special  computer 
programs,  such  as  GE  STATPAC,  see  Nelson,  Morgan  and 
Capcral  L11J  which  are  often  not  readily  available  to 
analysts.  Also,  maximum  likelihood  estimates  are  subject 
to  small  sample  bias.  Graphical  methods  provide  useful 
supplements  to  a more  formal  analysis,  but  are  subjective 
in  nature.  The  linear  estimation  methods  are  computationally 
simpler  than  the  maximum  likelihood  methods.  However, 
unlike  the  procedures  we  propose  in  this  paper,  they  are 
strictly  applicable  to  Type  II  censoring  only.  This  means, 
for  example,  that  the  information  at  test  conditions  where 
no  failures  have  occurred  is  ignored  in  the  analysis.  These 
conditions,  however , are  frequently  the  ones  of  greatest 
practical  interest. 

Hartley  t 4 ) proposes  an  iterative  maximum  likelihood 
approach  which  bears  some  similarity  to  the  proposed 
iterative  least  squares  procedure.  In  Hartley's  paper, 
the  missing  observations  are  replaced  by  introducing 
pseudo-frequencies  and  applying  the  standard  maximum 
likelihood  procedures  for  complete  data. 
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Further  discussion  of  this  method  and  some  extensions 
are  given  by  Hartley  and  Hocking  I.  5 1,  Hughes  f 6 1,  and 
Krane  T 8 ].  In  a recent  paper,  Dempster,  Laird  and 
Rubin  II  2 1 describe  an  algorithm  for  maximum  likelihood 
estimation;  this  algorithm  can  also  be  applied  to  the 
problem  considered  in  this  paper  (see  Section  7). 


3.  DESCRIPTION  OF  METHOD  FOR  SIMPLE  SITUATION 

Assume  the  standard  simple  regression  model  within 
the  region  of  interest  between  the  stress  x and  the 
average  time  to  failure  y for  some  device  on  life  test 
at  that  stress,  i.e., 

^x  = 30+BlX  (1) 

where,  at  any  stress,  (log)  time  to  failure  is  normally 

distributed  with  a constant  standard  deviation  a and  where 

Bq,  (3-^  and  o are  unknown  parameters.  One  or  more  units 

are  tested  at  each  of  several  stresses.  The  resulting 

data,  after  c^  units  of  running  time  at  the  various  stresses, 

consist  of  the  failure  times  (<c  ) on  the  failed  units  and 

— x 

the  running  times  c on  the  unfailed  units  or  "run-outs." 

x 

At  some  (usually  high  stress)  conditions,  there  may  be 
only  failures  and,  at  other  (usually  low  stress)  conditions, 
there  may  be  only  run-outs. 

Let  c denote  the  censoring  time  for  a particular 

I i 

run-out  at  stress  x.  The:n,  using  the  well-known  properties 
of  the  truncated  normal  distribution  (see  Johnson  and 
Leone  l 7 .]  or  Nelson  T9]),  the  expected  value  y*  of  the 
failure  time  for  this  unit  is: 

y*  = y + of  ( z)  / [l-F  ( z)]  (2a) 

X X 

where 

z = (c  - y )/o  (2b) 

X X 

and  f(z)  and  F(z)  denote,  respectively,  the  ordinate  at  z 
and  the  area  to  the  left  of  z of  a unit  normal  distribution. 
Thus,  for  this  situation,  the  iterative  least  squ ires 

I I 

procedure  is  as  follows: 
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Iteration  0 


Step  1:  Fit  the  linear  relationship  using  standard  least 

squares  regression  analysis,  treating  the  run-outs  as  if 

~ (0)  ~ (0) 

they  failed  a,_  their  censoring  times.  Let  Bq  , B-^ 
and  c(0)  denote  the  resulting  initial  estimates  of 
Bq,  B1  and  o respectively. 

Step  2:  Use  the  initial  regression  fit  from  Step  1 to 
estimate  the  unconditional  mean  time  to  failure  for  each 

run-out.  Denote  the  resulting  estimate  at  stress  x as 

"'(O)  ~(0)  ^ (0) 

yx  . Then,  use  the  initial  estimates  y^  and  a , 

instead  of  the  unknown  values  y and  a,  to  estimate  the 

x ' 

"*  ( 0) 

mean  failure  time  y^  for  each  run-out  x using  equations 
( 2a)  and  ( 2b) . 

Iteration  1 

Step  1:  Obtain  a revised  least  squares  regression  fit 

(o) 

using  the  estimated  mean  failure  times  y^  obtained  in 
Step  2 of  the  previous  iteration  for  the  run-outs.  Let 

~ (i)  ~ (i)  ~ (i) 

Bq  , 8-^  and  o denote  the  new  estimates  of  8q,  B-^  and  o. 
Step  2:  Repeat  Step  2 of  the  previous  iteration  using  the 
least  squares  estimates  obtained  in  Step  1 of  the  current  itera- 
tion to  re-estimate  the  mean  times  to  failure  for  the  run-outs. 
Subsequent  Iterations:  Continue  the  above  procedure  until 
convergence  is  achieved. 

As  indicated,  the  preceding  iterative  least  squares 
procedure  requires  only  a standard  computer  program  for 
regression  analysis  end  tabulations  or  calculations  of  the 
ordinates  and  tail  areas  for  a standard  normal  distribution 


or  a tabulation  of  the  expected  values  of  a truncated 
normal  variate.  In  practice,  one  might  wish  to  mechanize 
the  process  by  writing  a computer  routine  to  perform  the 
iterations.  We  have  developed  such  a program,  using 
standard  approximations  for  normal  distribution  ordinates 
and  areas.  This  program  is  used  for  the  simulation  analysis 
described  in  Section  5.  It  is  less  complicated  than  the 
companion  program  used  in  Section  5 for  obtaining  estimates 
by  the  method  of  maximum  likelihood. 

The  speed  of  convergence  depends  on  the  proportion 
of  censored  observation.  The  iterations  in  the  simulation 
analysis  were  stopped  after  both  the  slope  and  intercept 
estimates  agreed  to  3 decimal  places  on  two  consecutive 
iterations.  In  most  cases  convergence  was  achieved  in 
fewer  than  15  iterations  and  in  very  few  cases  were  more 
than  50  iterations  required. 
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4.  NUMERICAL  EXAMPLE 


Table  1 gives  the  results  of  temperature  accelerated 
life  tests  on  electrical  insulation  in  40  motorettes, 
originally  reported  by  Crawford  Cl  3.  Ten  motorettes  were 
tested  at  each  of  four  temperatures.  Testing  was  termi- 
nated at  different  times  at  each  temperature,  resulting  in 
a total  of  17  failed  units  and  23  unfailed  ones.  The  model 
used  to  analyze  the  data  assumes  that: 

i)  for  any  temperature,  the  distribution  of  time 
to  failure  is  lognormal 

ii )  the  standard  deviation  a of  the  lognormal  time 
to  failure  distribution  is  constant,  and 

iii)  the  mean  of  the  logarithm  of  the  time  to  failure 
Mx  is  a linear  function  of  the  reciprocal 
x = 1000/ (T+273 . 2)  of  the  absolute  temperature 
T,  that  is 

px  = 30  + exx 

where  Bq,  B-^  and  o are  unknown  parameters.  The 
preceding  model  is  often  referred  to  as  the 
Arrhenius  relationship. 

Nelson  and  Hahn  C 1 0 3 fitted  the  given  data  to  the 
preceding  relationship  using  linear  estimation  methods, 
and  Hahn  and  Nelson  T 3 3 used  graphical  and  maximum 
likelihood  estimation  techniques  to  fit  the  same  data 
and  compared  the  results  from  the  various  methods.  We 
now  fit  the  data  using  the  iterative  least  squares  method. 

The  procedure  is  described  below  and  illustrated  in  Figure  1. 
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Iteration  0 


Step  1:  The  censored  observations  are  initially  taken 
as  failures,  i.e.,  the  10  unfailed  units  at  150°C  are 
assumed  to  be  failures  at  8064  hours,  the  3 unfailed  units 
at  170°C  are  assumed  to  be  failures  at  5448  hours,  etc.  Then 

simple  least  squares  is  used  to  fit  the  assumed  relation- 

ship. This  leads  to  the  following  initial  parameter 
estimates : 

B0(0)  = -4.9307,  81(°)  = 3.7471, 
o(0)  = 0.1572 

Step  2:  From  the  initial  fit,  unconditional  expected 
failure  times  for  each  of  the  censored  values  are  estimated. 
For  example,  at  220°C,  i.e.,  (x  = 1000/(220  + 273.2)  = 2.03), 

the  mean  log  time  to  failure  is  estimated  from  the  initial 

fit  to  be 

W(0)  2.03  = -4.9307  + 3 . 7 47 1 [l000/ ( 220+273 . 2)]  = 2.6668. 

Thus,  based  upon  the  initial  fit,  log  time  to  failure  at 
220°C  is  estimated  to  follow  a normal  distribution  with  an 
estimated  mean  of  2.6668  and  an  estimated  standard  deviation 
of  0.1572. 

Censoring  at  220°C  occurred  at  528  hours,  i.e., 
c 2 03  = log  (528)  = 2.7226 

If  one  assumes  a log  normal  distribution  for  time  to 
failure,  the  censoring  time  is  equivalent  to  an  estimated 
normal  deviate  of 


z = 2.03 


2.03  = 2.7226  - 2.6668  * 0.3552 
0.1572 
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The  resulting  estimate  of  the  expected  log  time  to  failure 
for  each  of  the  5 unfailed  units  at  220°,  conditional 
upon  their  observed  running  times  of  528  lours,  is 
therefore 


(0) 


(0) 


(0) 


*J  2 . 0 3 M2.03 


+ 0 I^FBT  ■ 2-6668  + <0-15721  1 -h'.iltl  ‘ 2 

or  676  hours.  Similarly,  estimates  of  the  conditional 
expected  failure  times  for  the  unfailed  units  at  the  other 
tlree  temperatures  are 

150°C:  4.0384  log  hours  (or  10,924  hours) 

170°C : 3.8089  log  hours  (or  6,440  hours) 


0.3745 


190°C:  3.3295  log  hours  (or  2,135  hours) 


Iteration  1 

Step  1:  The  estimated  log  times  to  failure  for  the 
censored  observations  obtained  in  Step  2 of  Iteration  0 
are  used  to  fit  a new  least  souares  line  to  the  data. 

The  following  parameter  estimates  result: 

B0(1)  = -5.2603,  31(1)  = 3.9263,  o(1)  = 0.1799 
Step  2:  Using  the  estimates  from  Step  1 of  this  iteration 
the  following  revised  estimate  of  the  conditional  expected 
time  to  failure  is  calculated  for  each  of  the  censored 
observations : 

150°C : 4.09852  log  hours  (or  12,546  hours) 

170°C : 3.83972  log  hours  (or  6,914  hours) 

190°C : 3.36553  log  hours  (or  2,320  hours) 

220  C:  2.85846  log  hours  (or  722  hours)  =11^  Q3 
Subsequent  Iterations:  The  preceding  steps  are  repeated 
until  convergence  is  achieved.  This  example  required  17 


.83 


12 


1 


iterations.  The  final  parameter  estimates  are 

A A A 

B0  = - 5.81829,  = 4.20426,  o = 0.204322 

and  the  final  estimates  of  the  expected  time  to  failure  for 
the  censored  observations  are 

150°C : 4.17297  log  hours  (14,893  hours) 

170°C:  3.87676  log  hours  (7,529  hours) 

190°C : 3.40090  log  hours  (2,517  hours) 

n * ( 17) 

220  C:  2.87982  log  hours  (758  hours)  = u 2 
Table  2 shows  the  results  obtained  at  various 

iterations.  In  this  example  results  close  to  the  final 

values,  from  a practical  viewpoint,  were  obtained  after 

12  iterations. 

The  iterative  least  squares  estimates  are  compared 
in  Table  3 with  the  maximum  likelihood  estimates  obtained 
by  Hahn  and  Nelson  C33  for  the  same  data.  The  results 
from  the  two  methods  are  quite  similar  and,  except  for 
the  residual  error  standard  deviation,  agree  more  closely 
with  each  other  than  with  the  estimates  obtained  from 
graphical  and  linear  unbiased  estimation.  The  latter 
two  methods  do  not,  however,  use  the  test  data  at  150°C, 
where  none  of  the  10  test  units  failed.  If  one  ignores 
the  150°C  data,  the  iterative  least  squares  estimates 
and  the  maximum  likelihood,  estimates  were  even  closer  to 
each  other  than  before  and  again  agreed  better  with  each 
other  than  they  did  with  the  graphical  and  linear  unbiased 
estimates  (details  available  from  the  authors) . 
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RESULTS  OF  SIMULATION  ANALYSIS 
A . Description  of  Study 

A simulation  analysis  was  used  to  compare  the 
properties  of  the  iterative  least  squares  procedure 
with  those  of  the  method  of  maximum  likelihood.  A 
linear  relationship  was  assumed  between  the  mean  (log) 
time  to  failure  and  the  applied  (log)  stress  x , i.e., 

ux  = 30  + 31X 

and  at  a given  stress,  the  (log)  time  to  failure  was 
assumed  to  follow  a normal  distribution  with  a constant 
standard  deviation  a , where  8Q,  8-^  and  o are  unknown 
parameters . 

Without  any  loss  of  generality,  8q  and  8^  were  set 
equal  to  1 and  -1,  respectively.  Testing  was  assumed 
to  be  conducted  uniformly  over  the  interval  (0,1),  i.e., 
ranging  from  the  minimum  stress  of  x = 0 (with  y^  = 1) 
to  the  maximum  stress  of  x = 1 (with  y^  = 0). 

Simulation  analyses  were  conducted  at  each  of  the 
14  different  conditions  shown  in  Table  4.  The  following 
variables  were  studied: 

1.  The  censoring  probability  p at  the  minimum 
stress,  i.e.,  at  x = 0. 

2.  The  ratio  r of  the  censoring  probability  at 
the  minimum  stress  (x  = 0)  to  the  censoring 
probability  at  the  maximum  stress  (x  = 1) . 

3.  The  number  of  equally  spaced  test  stresses  k 
over  the  interval  from  x = 0 to  x = 1 . One  test 
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unit  per  stress  was  assumed  at  11  of  the 
conditions;  multiple  test  units  per  stress 
were  evaluated  at  the  other  three  conditions. 

Specification  of  p and  r,  together  with  the  stated 
assumptions,  defines  the  censoring  time  c^  at  each 
stress  and  also  the  common  stanoard  deviation  o . 

The  14  conditions  were  as  follows: 

- The  "center"  condition  of  p = 0.75,  r = 5 and 

k = 10;  thus,  the  censoring  probabilities  over  the 
10  stresses  ranged  from  0.75  at  x = 0 to  0.15  at 
x = 1 (Condition  1) . 

- The  8 combinations  of  conditions  (in  a full  factorial 
arrangement)  of  p = 0.5  and  0.95,  r = 2 and  10  and 

k = 5 and  20  (Conditions  2 through  9) . 

- The  two  conditions  with  p = 0.75  and  r = 5 (in  both 
cases)  and  with  k = 5 and  20  (Conditions  10  and  11) . 

- The  condition  with  p = 0.75,  r = 5 and  3 test  units 
at  x = 0 and  x = 1 and  4 test  units  at  x = 0.5 
(Condition  12) . 

- The  condition  with  p = 0.75,  r = 5 and  5 test  units 
at  x = 0,  3 test  units  at  x = 0.5  and  2 test  units 
at  x = 1 (Condition  13) . 

- The  condition  with  p = 0.75,  r = 5 and  2 test  units 
at  x = 0,  3 test  units  at  x = 0.5  and  5 test  units 
at  x - 1 (Condition  14) . 

Enough  simulations  were  conducted  at  each  of  these 
14  conditions  so  as  to  obtain  1,000  runs  (5,000  and 
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10,000  runs  at  some  conditions)  with  one  or  more 
censored  observations  and  two  or  more  uncensored 
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ones.  Table  4 shows  for  each  condition  the  number 
of  evaluated  simulation  runs,  the  average  number  of 
censored  observations  in  these  runs,  and  also  percent 
of  added  unevaluated  runs  (relative  to  the  number  of 
evaluated  runs)  at  each  test  condition  resulting  in 

0,  k-1  and  k censored  observations. 

In  each  simulation  run,  failure  times  were  randomly 
generated  under  the  previously  stated  assumed  model. 

The  assumed  model  was  fitted  using  least  squares 
regression  analysis  based  on  the  complete  data  (CLS) 
before  censoring.  All  observations  exceeding  c„  were 

X 

then  assumed  to  be  censored  at  c^  and  new  fits  were 
obtained  using  iterative  least  squares  (ILS)  and 
maximum  likelihood  (ML) . 

B . Summary  of  Results 

1 . Overall  Summaries 

Tables  5 and  6 summarize  the  following  results 
of  the  simulation  analyses: 

- The  estimated  regression  intercept  and  mean 

at  the  low  stress  (Bq  = Mq  = D (Tables  5A  and  6A) . 

- The  estimated  regression  slope  (B-^  = -1)  (Tables 
5B  and  6B) . 

- The  estimated  mean  at  the  high  stress 

(B0  + Bx  = = 0)  (Table  5C) . 

- The  estimated  standardized  residual  error 
standard  deviation,  i.e.,  the  ratio  of  the 
calculated  residual  error  standard  deviation 
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to  the  true  residual  error  standard 
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deviation  (o'  = 1)  (Tables  5D  and  6C) . 

Each  part  of  Table  5 provides  a comparison, 
for  each  test  condition,  of  (1)  the  average, 

(2)  the  standard  deviation,  (3)  the  root  mean 
square  error  and  (4)  the  minimum  and  maximum 
values  of  the  desired  estimates  using  the  CLS , 

ILS  and  ML  methods.  This  table  is  based  upon 
those  simulation  runs  which  resulted  in  one  or 
more  censored  observations  and  two  or  more 
uncenscred  ones.  The  simulations  which  did  not 
result  in  any  censored  observations  were  excluded 
because,  in  these  cases,  CLS,  ILS  and  ML  provide 
identical  results  (except  that  the  ML  estimate 
of  the  residual  error  standard  deviation  is 
smaller  than  the  CLS  and  ILS  estimates  by  a 
factor  of  [(n-2)  /n]  Also,  those  simulations 

with  all  censored  observations  or  with  only 
a single  uncensored  observation  were  excluded 
because  no  ML  estimates  can  be  obtained  for  them. 
As  a result  of  these  exclusions,  the  resulting 
(conditional)  CLS  estimates  of  Bg  and  are 
no  longer  unbiased. 

Table  6 provides,  for  each  test  condition, 

(1)  paired  comparisons  of  the  average  difference, 

(2)  standard  deviation  of  differences  (3)  mini- 
mum and  maximum  differences  in  estimates  for 
ILS  - ML,  ILS  - CLS,  and  ML  - CLS  with  regard 
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to  and  Rj  and  similar  statistics  for  the 
ratios  with  regard  to  o ' and 

(4)  evaluations  of  Pitman  closeness  (see  Pitman 
f 1 3 1 ) , i.e.,  the  percentage  of  simulation 
runs  for  which: 

- The  ILS  estimates  are  closer  to  the  true 
parameter  value  than  are  the  ML  estimates. 

- The  ILS  estimates  are  closer  to  the  true 
parameter  value  than  are  the  CLS  estimates. 

- The  ML  estimates  are  closer  to  the  true 
parameter  value  than  are  the  CLS  estimates. 

Specialized  Summaries 

Table  7 provides  detailed  two-way  tabulations 
showing  the  estimates  of  the  regression  intercept 
and,  equivalently,  the  mean  at  the  low  stress  (f?Q  = 1) 
from  the  5000  simulations  conducted  at  Condition  1 
for  the  following: 


- Iterative  least  squares  versus  maximum 
likelihood  (Table  7A) . 


- Iterative  least  squares  versus  complete 
least  squares  (before  censoring)  (Table  7B) . 

- Maximum  likelihood  versus  complete  least 
squares  (before  censoring)  (Table  7C) . 

Although  not  given,  similar  detailed  comparisons 

have  been  obtained  for  the  estimates  of  the  slope 
coefficient  and  of  the  residual  error  standard 
deviation,  and  also  for  each  of  the  other  13 


conditions . 
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c. 


Discussion  of  Results 


Inspection  of  Tables  5,  6 and  7 indicates  that, 
for  most  of  the  conditions  considered,  the  iterative 
least  squares  method  and  the  method  of  maximum 
likelihood  yield  very  similar  results  (see  especially 
Table  7A) . A notable  exception  was  Condition  6 
(discussed  further  below) . 

As  one  might  expect,  in  some  cases  both  ILS  and  ML 
gave  biased  estimates  of  appreciable  magnitude  at 
some  or  all  of  the  5 conditions  where  there  was  a 
large  number  of  unevaluated  simulations  due  to  obtaining 
0,  k-1  or  k censored  observations  (Conditions  2, 3, 4, 5 
and  10) . In  these  cases,  the  evaluated  CLS  estimates 
are  also  biased,  although  sometimes  the  magnitude 
of  the  bias  is  smaller.  In  addition  to  these  5 condi- 
tions, estimated  biases  exceeding  10%  resulted  in 

- Estimating  Bq  and  B^  using  both  ILS  and  ML  at 
two  further  conditions. 

- Estimating  Bq  + B1  (i.e.,  the  mean  at  the  high  stress 
condition)  using  ILS  at  one  further  condition. 

- Estimating  o'  using  ILS  at  seven  further  conditions 
and  using  ML  at  four  further  test  conditions. 

The  estimated  biases  with  regard  to  Bq,  B-^  and 
Bq  + B-^  were  positive  at  some  conditions  and  negative 
at  others.  In  evaluating  a' , however,  underestimates 
were  obtained,  on  the  average,  at  almost  all  condi- 
tions for  both  ILS  and  ML. 


It 
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An  overall  comparison  of  the  performance  of  the 
iterative  least  squares  and  maximum  likelihood 
estimation  methods  at  the  14  test  conditions  in 
the  simulation  analysis  is  shown  in  Table  8.  This 
tabulation  shows  that  there  were 

- An  approximately  equal  number  cf  conditions  where 
the  average  ILS  estimate  is  closer  to  the  true 
parameter  value*  than  is  the  average  ML  estimate. 

- A larger  number  of  conditions  at  which  the  root 
mean  square  eiror  is  smaller  for  the  ILS  estimate 
than  for  the  ML  estimate,  especially  with  regard 
to  Bq  and  B^. 

- A larger  number  of  conditions  at  which  the  ILS 
estimate  has  more  simulation  runs  closer  to  the  true 
parameter  value  than  does  the  ML  estimate  with 
regard  to  BQ,  B^  and  o'.  In  estimating  o',  however, 
the  average  percentage  of  runs  (averaged  over  the 

14  conditions,  weighting  each  condition  equally) 
at  which  the  ILS  estimate  is  closer  to  o'  was 
essentially  equal  to  that  for  which  the  ML  estimate 
was  closer. 

On  balance,  then,  the  proposed  iterative  least  squares 
procedure  gave  results  which  are  at  least  competi tive 
with  those  obtained  by  maximum  likelihood.  Iterative 
least  squares  did  poorest  relative  to  maximum  likelihood 
at  Condition  6,  although  maximum  likelihood  results 
at  that  condition  were  also  quite  poor.  The  censoring 
probability  at  this  condition  ranges  from  0.95  at 
the  low  stress  to  0.475  at  the  high  stress  and,  on 
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the  average,  about  three-quarters  of  the  20  observa- 
tions are  censored,  the  highest  proportion  among  the 
14  test  conditions.  The  large  number  of  run-outs, 
especially  at  or  near  the  low  stress,  d)  greatly  bias 
the  initial  fit  and  this  is  not  fully  compensated 
for  in  the  later  iterations  and  (2)  lead  to  an 
appreciable  underestimate  of  the  residual  error 
standard . 

An  underestimate  of  the  residual  error  standard 
deviation  estimated  by  the  iterative  least  squares 
method  versus  that  from  maximum  likelihood  also  was 
found  at  two  other  conditions  (Conditions  8 and  14) 
at  which  there  are  multiple  censored  observations  at 
the  same  or  neighboring  stresses.  (This  might  also 
have  been  the  case  in  the  motorette  example  of 
Section  3.)  Since,  the  censored  observations  at  the 
same  (or  similar)  stresses  have  the  same  (or  similar) 
estimates  of  the  expected  value,  the  variability  among 
them  is  small  and  thus  contributes  little  to  the 
overall  estimate  of  variability.  Some  possible  remedies 
for  this  are  suggested  in  Section  7. 

In  conclusion,  we  note  that  the  relatively  good 
results  of  iterative  least  squares  when  compared  with 
maximum  likelihood  may  be  due  to  the  small  sample  bias 
in  the  maximum  likelihood  estimates.  An  assessment  of 
the  magnitude  of  this  bias  is  a by-product  of  our 
simulation  analysis. 
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D.  Comment  on  Calculational  Requirements 


A major  reason  for  proposing  the  iterative  least 
squares  approach  is  that  the  calculations  are  less 
than  ttose  required  for  maximum  likelihood  estimation. 

For  example,  the  iterative  least  squares  calculations 
can  be  programmed  with  fewer  than  5<J  FORTRAN  statements, 
whereas  the  maximum  likelihood  calculations  using  a 
Powell  optimization  routine  require  approximately  300 
FORTRAN  statements.  The  average  number  of  function 
evaluations  for  the  14  test  conditions  in  the  simulation 
analysis  was  about  15.  For  problems  of  this  type, 
maximum  likelihood  usually  requires  at  least  50  function 
evaluations.  It  is  possible,  however,  that  some  of 
the  more  recently  suggested  maximum  likelihood  algorithms 
(see,  for  example,  Dempster,  Laird  and  Rubin  [ 21)  reduce 
the  required  calculational  effort  for  this  method. 
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6. 


GENERALIZATION  OF  RESULTS 


Most  of  the  preceding  discussion  has  dealt  with  the 
special  situation  where: 

- There  is  a simple  linear  relationship  (between  stress 
and  time  to  failure) . 

- The  random  variation  follows  a (log)  normal  distribution. 

- Censoring  is  to  the  right,  i.e.,  the  censored  values 
equal  or  exceed  the  observed  values  and  occur  at  the 
same  value  (or  time)  for  a particular  condition  (or 
stress) . 

The  iterative  least  squares  procedure  is  not,  however, 
limited  to  this  specialized  situation. 

More  general  regression  relationships 

The  proposed  method  handles  censoring  of  the  dependent 
variable.  It  is  not  limited  to  any  special  form  of  relation- 
ship and  can  be  applied  to  multiple  regression  or,  more 
generally,  to  any  situation  for  which  least  squares  pro- 
cedures are  used. 

Other  distributions  for  random  variation 

Let  f(yx)  and  F (y  ) denote  the  probability  density 
function  and  cumulative  distribution  function,  respectively, 
of  the  random  variable  y at  the  condition  x and  let  censoring 
to  the  right  occur  at  the:  value  c^  of  y.  Also  let  y*  denote 
the  value  which  would  be  observed  for  a censored  observation 
if  it  were  observable,  e.g.,  the  time  at  which  a run-out 
would  fail  if  there  were  no  censoring.  Then  the  conditional 
probability  density  function  of  an  observation  censored  to 
the  right,  i.e.,  one  which  is  known  to  exceed  the  censoring 
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value  c , is  that  of  the  random  variable  y truncated  to 

X X 
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the  riqht  at  c . Therefore 
x 

f<rl>  - f‘yx)/[l-F(Cx)]  , yj  > cx 

and  the  conditional  expected  value  of  the  censored  observa- 
tion (i.e.  its  expected  failure  time)  is 

“I  = ;c  (yx)dyx  /h-FfSt)]  • (3) 

x 

Thus,  the  proposed  method  requires  evaluation  of  the  conditional 
expectation  given  in  equation  (3).  If  the  random  variation 
is  normally  distributed  this  results  in  equations  (2a)  and 
(2b).  Similar  results  can  be  obtained  for  ether  distributions 
for  the  random  variation. 

Other  censoring  schemes 

In  practice,  one  may  encounter  situations  involving 
censoring  other  than  to  the  right.  For  example,  in  a life 
test,  censoring  may  occur  at  random  times  due  to  accident  or 
breakdown  of  equipment.  Also,  cen  soring  other  than  to  the 
riqht  can  occur  when  failures  are  discovered  only  at  periodic 
inspections . 

The  iterative  least  squares  concept,  just  as  the  method 
of  maximum  likelihood,  applies  to  such  more  general  censoring 
schemes.  However,  when  the  censoring  occurs  early,  it 
might  be  better  in  the  initial  least  squares  fit  to  ignore 
the  censored  observations  instead  of  treating  them  as 
failures.  Also,  a procedure  which  assigns  varying  lesser 
weights  to  the  censored  observations  than  to  the  uncensored 
ones  (see  next  section)  might  be  especially  appropriate  for 
more  general  censoring  schemes. 
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7.  AREAS  FOR  FURTHER  STUDY 


A theoretically  based  study  of  the  properties  of  the 
proposed  iterative  least  squares  approach  is  warranted.  In 
addition,  some  specific  areas  for  further  consideration  are 
reviewed  below. 

Assigning  less  weight  to  censored  observations 

The  informational  content  of  the  censored  observations 
is  less  than  that  of  the  uncensored  ones.  For  example, 
in  censoring  to  the  right,  the  earlier  the  censoring  occurs, 
the  less  is  the  information  provided  by  the  censored  observa- 
tion. This  might  be  taken  into  consideration  in  the 
iterative  least  squares  analysis  by  assigning  less  weight 
to  the  censored  observations  than  to  the  uncensored  ones  and 
by  varying  these  weights  depending  on  the  censoring  point. 

The  following  is  one  possible  weighting  scheme  for  tie  censored 

observation  v*  at  the  condition  x: 
x 

w^x)  = 1 - [o2(v*)/o2] 

where  a denotes  the  (constant)  residual  error  standard  devia- 
tion and  o(y*)  denotes  the  standard  deviation  of  the  truncated 
x 

distribution  of  an  observation  known  to  exceed  c . For 
example,  in  a product  life  analysis  with  censoring  to  the 
right,  a denotes  the  standard  deviation  of  the  distribution 
of  time  to  failure  at  a constant  stress  and  o( y*)  denotes 
the  standard  deviation  of  the  distribution  of  the  time  to 
failure  of  a unit  at  stress  x which  is  known  tc  be  unfailed 
at  time  cy . In  particular,  for  normally  distributed  random 
variation,  it  can  be  shown  (see  Nelson  l9J)  that,  using  the 
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notation  of  Section  3, 

, ,1/2 

o ( y* ) = all  + [_zf  (z)  / f 1 - F ( z )]J  - [f(z)/(l  - F(z))Jz} 
where,  as  before,  z = (c  - y ) /a . 

X X 

This  procedure  assigns  a weight  of  1 to  all  uncensored 
observations  and,  for  normally  distributed  random  variation, 
leads  to  a weight  which  approaches  0 with  decreasing  c^  and 
a weight  of  1 with  increasing  c^.  This  scheme  may,  however, 
be  inappropriate  when  the  random  variation  is  not  normally 
distributed.  For  example,  if  the  time  to  failure  distribu- 
tion is  exponential,  all  the  censored  observations  would 
receive  a weight  of  0.  For  a Weibull  distribution  with  a 
shape  parameter  less  than  1,  censored  observations  would 
receive  negative  weights. 

An  alternative  simple  scheme  is  to  weight  the  censored 
observation  proportional  to  the  associated  cumulative 
(failure)  probability  at  the  censoring  point,  i.e.,  to 
assign  to  an  observation  censored  to  the  right  of  c^  the 
weight 

w2(x)  = F(cx) 

and  to  give  all  uncensored  observations  a weight  of  1. 

In  practice,  the  required  distribution  ordinates  and 
cumulative  probabilities  are,  of  course t unknown.  In 
their  place  one  would  use  the  estimates  obtained  in  the 
most  recent  iteration.  For  the  initial  iteration,  some 
working  rule  needs  to  be;  developed  to  obtain  an  appropriate 
(initial)  weight  between  0 and  1 for  the  censored  observations. 

The  reason  for  using  weighted  iterative  least  squares 
estimation  is  the  expectation  that  it  will  lead  to  better 


performance  than  the  proposed  unweighted  estimator.  The 
weighting  procedure,  however,  complicates  the  method  and 
this  compromises  our  claim  of  simplicity.  Despite  this, 
a computer  program  based  upon  a weighted  iterative  least 
squares  approach  is  still  likely  to  be  appreciably  simpler 
to  develop,  use  and  explain  than  one  based  upon  maximum 
likelihood . 

Development  of  statistical  inference  procedures  and 
modification  of  standard  regression  programs 

The  preceding  discussion  has  been  concerned  with 
obtaining  point  estimates  using  the  proposed  iterative  least 
squares  procedure.  Frequently,  statistical  intervals  and 
hypothesis  tests  are  also  desired.  With  the  method  of 
maximum  likelihood,  this  is  accomplished  using  asymptotic 
theory.  Such  inferences  tend  to  be  unconservative  for 
small  sample  size  situations. 

Using  an  iterative  least  squares  approach,  one  can 
obtain  unconservative  inferences  by  applying  the  usual 
least  squares  inference  procedures  to  the  results  of  the 
final  iteration  and  treating  the  censored  observations  as 
if  they  were  uncensored.  In  addition,  it  is  possible  to 
obtain  conservative  confidence  intervals  and  prediction 
intervals  by  ignoring  the:  censored 

observations  altogether  in  obtaining  the  length  of  the 
desired  interval  (but  not  in  obtaining  the  point  estimate 
around  which  the  interval  is  constructed)  after,  perhaps, 
adjusting  for  bias  in  the  estimate  of  the  residual  error 


standard  deviation. 
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A compromise  between  such  overly  unconservative  and 
overly  conservative  approximations  would  be  to  calculate 
the  bounds,  as  well  as  the  original  estimates,  using  a 
scheme  which  assigns  varying  v;eights  to  the  censored  observa- 
tions (see  above  discussion). 

Along  similar  lines,  consideration  might  be  given 
to  appropriate  ways  of  adapting  standard  computer  programs 
for  regression  analysis  to  handle  residual  plots  with 
censored  data  and  for  performing  various  standard  analysis 
of  variance  tests  to  assess  the  significance  and  adequacy 
of  a fitted  model  with  censored  data.  In  particular,  a 
referee  has  suggested  a lack  of  fit  test  based  on  censored 
versus  non-censored  errors  as  a natural  extension  of  our 
procedures  and  as  a useful  method  of  evaluating  the  lack  of 
fit  of  the  assumed  equation  (see  also  Nelson  1 1 2l  ) . 

Specialization  to  a single  condition 

The  iterative  least  squares  procedure  is  not  limited 
to  situations  dealing  with  a regression  relationship.  It 
can  similarly,  and  more  simply,  be  applied  when  data  at  a 
single  condition  are  available  and  estimates  are  desired 
only  for  the  parameters  of  the  distribution  for  the  random 
variation  at  that  condition.  We  have  been  concerned  with 
the  more  general  regression  situation  because  (i)  this  is 
the  type  of  practical  problem  which  we  have  encountered 
most  frequently,  and  (ii)  special  simplified  procedures 
are  already  available  for  the  single  condition  situation. 

However,  this  situation  might  be:  useful  for  studying  various 
generalizations  and  extensions,  such  as  the  properties 
of  iterative  least  squares  under  different  censoring 
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schemes  or  the  effect  cf  different  start-up  or  weighting 
procedures.  The  simpler  single  condition  situation  might 
also  be  appropriate  for  a theoretical  study  of  the  properties 
of  the  proposed  procedures. 

Use  of  results  from  early  iterations 

The  number  of  iterations  to  convergence  is  important 
for  those  who  use  an  iterative  least  squares  approach  to 
analyze  censored  data  via  a standard  least  squares  computer 
routine  or,  for  example,  a hand  calculator.  It  is  also  of 
interest  to  determine  how  close  one  can  come  to  the  final 
convergence  values  after  a small  number,  such  as  3,4,  or 
5,  iterations.  In  this  regard,  procedures  to  speed  conver- 
gence also  merit  consideration. 

Bias  adjustment  procedures 

The  results  of  the  simulation  analysis  indicate  bias 
at  some  of  the  conditions  in  the  iterative  least  squares 
estimates,  as  well  as  in  the  maximum  likelihood  estimates, 
especially  in  the  estimate  of  the  residual  error  standard 
deviation.  The  simulation  results  might  themselves  be  used 
to  obtain  appropriate  bias  correction  factors. 

Procedures  for  correcting  for  the  underestimate  of  the 
residual  error  standard  deviation  within  the  iterative  scheme 
also  warrant  consideration.  An  upward  adjustment  factor 
might  be  arrived  at  based  upon  the  proportion  of  censored 
observations  in  the  data.  One  possibility  is  to  reduce  the 
degrees  of  freedom  in  calculating  the  residual  error  standard 
deviation.  In  a private  communication,  R.  Regal  has 
suggested  to  us  an  alternative  scheme  for  calculating  the 
residual  error  standard  deviation  based  upon  the 


29  - 


EM-algorithm  discussed  in  a paper  by  Dempster, 

Laird  and  Rubin  C 2 L The  EM  procedure  can  be  viewed  as 


a generalization  of  the  method  proposed  by  Hartley  C41  and 

by  other  authors  referenced  in  L21  . It  uses  the  same  expression  for 

the  slope  and  intercept  estimates  as  our  proposed  method. 


8.  CONCLUSION 


This  paper  proposes  an  alternative  scheme  for  analyzing 
censored  data  which  might  be  more  attractive  to  some  practi- 
tioners than  existing  methods.  This  is  because  the  proposed 
procedure  uses  well-known  least  squares  techniques,  requires 
less  computation  than  maximum  likelihood,  and  should  be 
easier  to  explain  to  non-statisticians. 

A simulation  analysis  indicates  that  the  proposed 
method  gives  statistically  reasonable  results  which,  for 
the  situations  considered  in  the  evaluations,  are  at  least 
competitive  with  those  from  maximum  likelihood.  Although 
the  simulations  were  limited  to  cases  where  censoring  is 
to  the  right  and  where  the  censored  values  exceed  the 
uncensored  ones,  this  is  the  situation  which  is  frequently 
encountered  in  the  analysis  of  accelerated  life  data  and 
many  other  practical  applications. 

As  indicated  in  Sections  6 and  7,  various  generaliza- 
tions of  these  results  are  possible.  To  extend  and 
refine  the  results  obtained  to  date,  further  work  is 
warranted . 
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APPLICATION  OF  ITERATIVE  LEAST  SQUARES  METHOD  TO 
INSULATION  LIFE  ACCELERATED  TEST  PROGRAM  DATA 


TABLE  1 


INSULATION  LIFE  IN  HOURS  AT  VARIOUS  TEST  TEMPERATURES 


150°C 


1) 


Test  Temperature 


170°C2) 

31 

190°C  ' 

220°C 

1764 

408 

408 

2772 

408 

408 

3444 

1344 

504 

3542 

1344 

504 

3780 

1440 

504 

4860 

5196 

^ All  10  motorettes  at  150°C  still  on  test  without 
failure  at  8064  hours. 

2^3  motorettes  at  170°C  still  on  test  without  failure 
at  5448  hours. 

^5  motorettes  at  190°C  still  on  test  without  failure 
at  1680  hours. 

^5  motorettes  at  220°C  still  on  test  without  failure 
at  528  hours. 


f 


TABLE 


TABLE  3 


COMPARISON  OF  RESULTS  FROM  ITERATIVE  LEAST  SQUARES 
AND  MAXIMUM  LIKELIHOOD  FITS  FOR  INSULATION  LIFE  DATA  EXAMPLE 


Estimates  from 


Iterative 

Least 

Maximum  Likelihood 

Squares 

Method 

Intercept  Coefficient  ($q) 

Slope  Coefficient  (g  ) 

-5.818 

-6.027 

A. 204 

(+1.489*) 

4.314 

1 

Log  Standard  Deviation  (o) 

0.204 

(+0.684*) 

0.259 

Median  Life  in  Hours  at: 

220°C 

508 

530 

190°C 

1812 

1940 

170°C 

4654 

5080 

150°C 

13060 

14680 

130°C*** 

40638 

47000 

10%  Point  of  Life  Distribution 
130°C*** 

at22230 

(27500  to  79500**) 

21900 

*Bounds  of  approximate  90%  confidence  interval 

**Approximate  90%  confidence  interval 

***Extrapolated  beyond  range  of  data 
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CLS  " Complete  least  squares  (before  censoring);  ILS  * Iterative  least  squares;  ML  - Maximum  likelihood. 
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Before  censoring 


TABLE  8 


OVERALL  COMPARISON  OF  THE  PERFORMANCE  OF  THE  ITERATIVE  LEAST  SQUARES 
AND  MAXIMIMUM  LIKELIHOOD  ESTIMATION  METHODS  AT  THE  14  TEST  CONDITIONS 

IN  THE  SIMULATION  ANALYSIS 


Number  of  Conditions  Out  of  14: 

Intercept 
and  Mean  at 
Low  Stress 
(B0) 

Slope 

(Bi> 

Mean 

at 

High 

Stress 

(Bo+bP 

Standardized 

Residual 

Error 

Standard 

Deviation 

(a’) 

1)  For  which  average  ILS  estimate 
is  closer  to  true  parameter 
value  than  i8  ML  estimate 

7 

6* 

4* 

7 

2)  For  which  ILS  estimate  has 
lower  root  mean  square  error 
than  ML  estimate 

12 

13 

9* 

8 

3)  For  which  ILS  estimate  has 
more  simulation  runs  closer 
to  true  parameter  value 
than  ML  estimate  has 

10 

10 

** 

10 

Average  percentave  of  runs  Cover 

14  conditions)  for  which  ILS 
estimate  is  closer  to  true 
parameter  value  than  ML  estimate 

* 

63.8% 

56.2% 

** 

50.6% 

Out  of  13  conditions. 
**Not  calculated. 
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times  are  estimated  for  the  censored  values.  These  are  then 
used  in  a further  least  squares  fit.  The  procedure  is  iterated 
until  convergence  is  achieved.  This  method  is  simpler  to  imple- 
ment and  explain  to  non-statisticians  than  maximum  likelihood 
and  appears  to  have  good  statistical  and  convergence  properties. 

The  method  is  illustrated  by  an  example,  and  some  simula- 
tion results  are  described.  Variations  and  areas  for  further 
study  also  are  discussed. 
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